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Abstract 

We present results of a numerical study of the differential equa- 
tion governing the stationary states of the two-dimensional planetary 
atmosphere and magnetized plasma (within the Charney Hasegawa 
Mima model). The most strinking result is that the equation appears 
to be able to reproduce the main features of the flow structure of a 
typhoon. 

1 Introduction 

There is a well known similarity between the two-dimensional models of the 
planetary atmosphere and the magnetized plasma. In the absence of dissi- 
pation the models can be reduced to differential equations having the same 
structure: the Charney equation for the nonlinear Rossby waves , in the 
physics of the atmosphere 0; and the Hasegawa-Mima equation for drift 
wave turbulence, in plasma physics |2j. They are similar with the Navier- 
Stokes equation because they have two conserved quantities, the energy and 
the enstrophy. This in principle allows states of negative temperatures, or, 
equivalently, these models support a trend to organised vortical flow. It re- 
sults the possibility to have as solutions coherent structures (vortices) besides 
the turbulent states characterised by spectral cascade. 
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These analytical models have led to a serious advancement of our knowl- 
edge in both fields. However the stationary states appear to be described 
within these models by a reduced equation having a too wide generality, 
representing actually something as a constraint with weak ability to iden- 
tify unequivocally the real solutions: it simply states that at stationarity the 
advection of the vorticity by the velocity vector field vanishes. In reality, 
numerical simulations show that the stationary states reached in relaxation 
are very regular and persist for a long time period and that this set of asymp- 
totic states is not the huge space of functions able to fulfill the constrained 
mentioned above. The fluid evolves at relaxation toward a reduced subset 
of functions, characterized by regular shape of the streamfunction PS]- 
|14j . [To"] (and references therein). At the oposite limit the turbulent regime 
can be treated with renormalization group methods [To] . 

It is well-known that the same phenomenon exists in the case of the ideal 
fluid described by the Euler equation. By experiments and numerical sim- 
ulation it has been shown that the ideal fluid evolves at relaxation toward 
a very ordered flow pattern, consisting of two (positive and negative) vor- 
tices and that this state persists for very long times, being limitted by only 
the effect of some residual dissipation. From numerical simulations it has 
also been inferred the form of the flow function. It has been found that 
the streamfunction obeys, in these states, the sm/i-Poisson equation. Mont- 
gomery and his collaborators have developed a theoretical statistical model 
which explains the appearence of this equation in this context [3], 0j, 0, 
M 0> 0) 0- Later, the equation has also been derived by formulating 
the continuum version of point-like vortices as a field theoretical model of 
interacting gauge and matter fields in the adjoint representation of SU (2) 
[17J. The essential point of the latter derivation was the self-duality of the 
relaxation states of the fluid. 

No equation (similar to the szn/i-Poisson equation in the Euler fluid case) 
has been found for the Charney-Hasegawa-Mima (CHM) equation, despite a 
considerable effort ^0], jTTj . However, as mentioned before, there are con- 
vincing experimental and numerical indications that the fluids (atmosphere 
and plasma) evolve to a reduced subset of states. 

We have developed a field theoretical model for the point-like vortices with 
short range interaction, based on Chern-Simons action for the gauge field in 
interaction with the nonlinear matter field, again in SU (2) algebra. It is 
then possible to derive the energy as a functional that becomes extremum on 
a subset of stationary states and presents particular properties. The general 
characterization of this family of states is their self-duality, which here means 
that the energy functional becomes minimum because the square terms are 
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all vanishing, leaving as lower bound a quantity with topological meaning. 
A very detailed account of the derivation is in Refs. [TH], [T§] . 

The result is a set of equations parametrized by the solutions of the 
Laplacean equation in two-dimensions. 

The simplest of these equations is 



(where p is a positive constant). There are already some confirmations that 
this is the equation governing the asymptotic stationary states of the CHM 
fluids : the scatterplots of (ip,cj) = (streamfunction, vorticity) obtained in 
experiments [20] and the scatterplots obtained in numerical simulations [10J 
are very similar to the nonlinear term of Eq. ((TJ . 

The objective of this work is to provide the first elements resulting from 
a numerical investigation of this equation. 

The results are summarised here. This differential equation is able to 
reproduce the main two-dimensional features of the typhoon vortical flow. In 
the physics of the atmosphere, it seems that other examples, like the tropical 
cyclones, can be reproduced by solutions of this equation. The following are 
the features we consider as very particular to the typhoon morphology (in 



1. The very narrow dip of the azimuthal velocity (mean tangential wind) 
in the center of the vortex, compared with the very large extension in 
space. This is characterized by the "radius of the maximum tangential 
wind" and this radius, as mentioned, is much smaller than the diameter 
of the vortex. Our equation is able to generate solutions with this 
structure. 

2. The slow decay of the magnitude of the azimuthal velocity toward the 
periphery, compared with the very fast decay toward the center; this is 
reproduced by the solutions of this equation. 

3. The very low magnitude (almost vanishing) of the vorticity over most 
of the vortex (approx. from the radius of maximum wind to the pe- 
riphery), while the magnitude in a narrow central region is extremely 
high. This feature is also reproduced by the equation. 

4. quantitatively, we obtain for the diameter of the typhoon's eye a rel- 
atively good magnitude. The vorticity is higher than in observations 
but not far from the realistic range. 



Aip H sinh ip (cosh ip — p) — 




2D) |2I], [22], [23, |2i: 
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We have very encouraging results of studies on plasma vortices, but they 
are not reported here. In plasma physics, the symmetrical, stable, vortical 
structures observed in experiments in the linear machine seem to belong to 
the class of solutions of this equation. We have also obtained several solutions 
that are very similar to the crystals of vortices, known from experiments. 

2 Numerical studies of the equation 

The numerical solution of this equation appears to be very difficult. This 
may be explained by the fact that the exponentials of the two functions sinh 
and cosh are very rapidly- varying functions and any perturbation is amplified 
and propagated in the solution. 

In addition, the Laplace operator has spurious solutions with exponential 
behavior that have to be eliminated by the numerical procedure. 

The paper of McDonald [2^] on the numerical integration of the sinh- 
Poisson equation is very helpful in understanding the problems related to 
a numerical treatment of our equation. However the approach proposed in 
that paper requires to use a small mesh, specifically for excluding the spuri- 
ous modes of the Laplacean. In the case of our equation, the vortices require 
a reasonable detailed description and this needs larger meshes. Then the 
problem of the precision of integration procedure arises and, if the initializa- 
tion happens to be far from one of the solution, the number of iteration of 
the solver is high and the errors accumulate, leading to lack of convergence. 
It may be supposed that the solutions would be similar to those of the sinh- 
Poisson equation, but structures with sharp spatial variation may be possible 

The structure of the function space representing the union of attractors 
for the various solutions of this equation appears to be very complex. This 
immediately translates into serious obstacles in the attempt to reach one of 
the presumed solution. The main instrument is, naturally, the initialization, 
i. e. to start the integration in the right subspace, representing the attractor of 
that solution. Since there is no available analytical description of this space, 
the search is simply a problem of guessing a reasonable initial function and 
to repeat as many times as necessary. One of the specific behaviors is the 
tendency of driving the solution toward the constant value 

i> = 4 l ' 2) (2) 

(see Eq.(j3J) which trivially verifies the equation. This seems to imply that 
there is a large attractor in the function space around these constant solu- 
tions. The solution which is larger in absolute magnitude is less stable since 
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any fluctuation around the constant generates high vorticity. We underline 
that the integrations described here are not radial (i.e. unidimensional) . 

With all the difficulties of getting a right initial positioning in the inte- 
gration procedure we note however that the solution with the typhoon mor- 
phology appears instistently from a wider class of initial shapes. 

2.1 The numerical code 

We use the code "GIANT A software package for the numerical solution of 
very large systems of highly nonlinear systems" written by U. Nowak and L. 
Weimann |2Zj. The code belongs to the numerical software library CodeLib 
of the Konrad Zuse Zentrum fur Informationstechnik Berlin. The 
meaning of the abbreviation is: GIANT = Global Inexact Affine Invariant 
Newton Techniques and corresponds to the implementation of the method 
proposed by Deuflhard (for many references see P7J). 
This code solves nonlinear problems 

F(x) = (3) 
initial guess of solution, x = xq 

The global affine invariant Newton schemes requires the solution of linear 
problems. For higher accuracy meshes the linear problems are solved by 
iterative methods. The balance between numerical requirements of the New- 
ton iteration (called outer iteration) and the iterative linear solver (inner) 
means that the solution of the linear problem will be approximative. Two 
packages of linear solvers can be used, GMRES (generalized minimum resid- 
ual : Brown, Hindmarsh, Seager) and GBIT1 (fast secant method using the 
Good-Broyden updates : Deuflhard, Freund and Walter). 

All necessary description of the method, of the code and many studies of 
the numerical precision and computer efficiency are presented by Nowak and 
Weimann in the documentation of the code. 

The code has been implemented and the tests have been performed with 
successful results (we are grateful to Dr. Weimann for his kind help in this 
problem) . 

2.2 Boundary conditions 

The boundary conditions are dependent on the value of p. The physical 
model imposes that the scalar function ip remains nonzero at infinity for 
p > 1. This means that we must require that the boundary condition is one 



5 



of the roots of the algebraic equation 



cosh-?/' — p = (4) 

which can give the vanishing of the physical vorticity at infinity. Then we 
impose 

boundary condition ip (r — > oo) = w b (5) 

= In (p ± \J p 2 — 

2.3 Initialization 

In general the initial profiles has been of two types: symmetric profiles with 
maximum centered on (0, 0) and initializations with functions expressed as 
product of trigonometric functions. 

The symmetric profiles has been chosen as Gaussian functions, or various 
annular shapes. 

For may runs, as suggested by the experiments for the sm/i-Poisson equa- 
tion (paper by McDonald the initial function is taken as a product of 
trigonometric functions in both directions, x and y. We need to prepare the 
initial function in the sense that the values that are obtained in for the vor- 
ticity, i.e. the Laplacean of the initial distribution should not be too different 
of what is obtained by simply inserting the initial function in the nonlinear 
term. For this we take a coefficient ipi n of the product of the trigonometric 
functions as a parameter to be determined. 

The initial function is taken as 

x/j (x, y) = ipjp + *l>in sin ( kir — — * mm j sin f kir— — ^ mm j (6) 

where k is the periodicity of the profile and ^ n is the amplitude. We insert 
in the equation and we require approximative equality of the two parts, the 
vorticity and the nonlinearity. This is obtained by choosing a point (x, y) 
where the initial function is maximum and it results a condition on only the 
amplitude, ip in . 

A^j = ^ m [2(k7t) 2 } (7) 
~ — sinh ip in (cosh ip in - p) 

This equation is solved and one of the roots is selected as the amplitude of 
the initial function. 
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The experiments with simple sin functions frequently lead to difficulties of 
convergence. Looking at the function's form (either partial evolutions during 
iterations or good, converged, results) we notice that the two-signed values 
are less tolerated and only one of the signs survives. This led us to adopt 
forms expressed as square of the trigonometric functions. 

3 Results of the numerical integration 

3.1 The typhoon morphology 

The value of the parameter is p = 1 . The domain is 

(x,y) E [-0.5,0.5] x [-0.5,0.5] 
with [101, 101] mesh points. The boundary value is 

tpP =ki(p- VV - l) = 
and the initial function is 

V (x, y) = 4 1] + 4.2 x sin (in % ~^ min ) sin (att V ~^ min 

It takes 501 calls to the function and Jacobian. The accuracy is 0.257 x 
10 -3 . This run has been executed with several mesh dimensions: [31 x 31], 
[51 x 51], [71 x 71]. The results are very close, but higher accuracy shows 
much clearer the details. 

The results are shown. The Figure ® shows the choice of the amplitude 
of the initialization and Fig.® shows the initial function if>. 

The solution has an apparent cylindrical symmetry around the center and 
for this reason we present a section along x of the streamfunction ip (x, y) 
(Fig.®). A section along x axis of the vorticity ui(x,y) is presented in 
Fig.®. 

In order to quantify the accuracy of integration we collect in all the do- 
main (x, y) the pairs (ip, to) and represent them together with the line of the 
nonlinear term in our equation, Fig.®. In Fig.® we show the ratio of the 
two quantities the nonlinear term and u, as resulted from the calculated ip. 
This ratio should be 1. There are points close to the value where this ratio 
is not 1 but, if we normalize adding an arbitrary constant to remove the 
possible singular cases, we notice a very good clustering of the points around 
the line 1. In addition, we represent the scatterplot of the pairs (u, magni- 
tudes of nonlinear term for the ^'s) and notice the close clustering around 
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Figure 1: The procedure to find an approximation to a good initialization. 





Figure 4: The vorticity, calculated from ijj(x,y) obtained by integration. 



the diagonal. Other tests are possible and they indicates that the integration 
is very good on most of the region and good within the imposed accuracy in 
the regions where the quantities reach values close to 0. 



Scatterptat (ip.W) ' rom integration. Tha line is the analytic nonlinear term 




0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 



Figure 5: The scatter plot (xjj, u), for p = 1. 



The contour plot of the solution is shown in Fig.fjHJ) on the same graph 
with the velocity field (we have used a reduced set of data due to limitations 
on the EPS file). We must note that this two-dimensional integration gives a 
radial component of the velocity which at maximum is about 20 times lower 
than the tangential one. 

The tangential component of the velocity is shown in two figures © and 
(jll)) with the purpose of making easier the observation of the central region. 
The narrow dip in the center is clearly visible and its radial extension can be 
compared with the extension of the whole domain. 

We have represented in Fig. (|T2*j) a section along the jC clXIS of the amplitude 
of the azimuthal component of the velocity. 
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Figure 6: The ratio of u and the nonlinear term. 



Scatlerplot (io,rh) compared with the diagonal 
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Figure 7: Scatterplot (u, the nonlinear term), compared with the diagonal. 
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Contours of streamfunction y(x,y) and velocity vector field (v ,v ) 
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Figure 8: The contours of the scalar streamfunction ip(x,y) and the vector 
field (v x .v y ). 
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The tangential component of the velocity, v , center is (0,0) 



20 




Figure 9: The tangential component vg(x,y) of the velocity vector field 
(v x .v y ), with center at (0,0). 
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The tangential component of the velocity, v center is (0,0) 





Figure 10: The tangential component vg(x,y) of the velocity vector field 
(v x .v y ), with center at (0,0) (same asEJ). 
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The tangential component of the velocity, v center is (0,0) 




Figure 11: The tangential component vg(x,y) of the velocity vector field 
(v x .v y ), with center at (0,0) (same asEJ). 
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The magnitude of the tangential velocity along a radius, v (r) 
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r 

Figure 12: The magnitude of the tangential component ve(x,y), seen along 
a radial line. The central fast decay is clearly visible. 
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3.1.1 Episodic structure of two vortices 

It is worth to mention that in a numerical experiment we have identified a 
state where two vortices have been formed, placed in symmetrical positions 
along the diagonal of the square domain [—0.5, 0.5] x [—0.5, 0.5] with a mesh 
of [31,31]. The value of the parameter is p = 1. The initial function is 
trigonometric with k = 2 in Eq.(0) with a coefficient ipu n = 3.8. It takes 
longer to obtain the solution with 0.84 x 10 -4 accuracy, 389 calls to the 
function. The result is in Fig.(|13)l. 

Contours of streamf unction iy(x,y) and velocity vector field (v^v J 




5 10 15 20 25 30 



Figure 13: The contours of the scalar streamfunction ip(x,y) and the vector 
field (v x .v y for a two- vortices approximative solution. 

This state has been reexamined with much higher accuracy. It has taken 
long time to see that the final solution was again the centered vortex shown 
before. Therefore from the point of view of the numerical experience this 
state of two vortices is irrelevant. However, the persistence of this state 
inside the iterative search may indicate that it is close to a solution, possible 
less stable. We have not investigated this further. Instead we will show below 
a solution with four vortices. 

3.1.2 Four vortices 

The calculations are done for p = 1 on meshes with various levels of details: 
31, 61, 101. The initial function is trigonometric with k = 3. 

The results show clearly the formation of four vortices, as shown by 
Fig. (fT5^1 . Each of them has a structure that is similar to the one presented 
in Fig.(jnj). It is interesting to note that again the vorticity is almost zero 
everywhere on the domain, except a strict region around the four vortices, 
where it reaches very high values. 
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To the accuracy we have used unitl now we cannot say if the local tan- 
gential velocity presents the same very fast decay to the center of the vortex. 

The solution streamfunction ^(x.y) 




Figure 14: The scalar streamfunction ijj(x,y) for a four-vortices solution. 

3.1.3 Four vortices obtained at p > 1 

For p = 3 it is also possible to obtain the four-vortex solution. The initial 
function is here a trigonometric combination for, k = 2 and squared such 
that only positive (four) maxima are initially present, with an amplitude of 
about ■0 = 4. 

3.1.4 The central strong decay of the tangential velocity, at p > 1 

The numerical integration is done for p = 3 , using an initialization by a 
centered peak from an trigonometric function. 

We note from Fig. (fTo^) that for larger values of the parameter p there is 
a even more narrow zone where there is the strong decay of the tangential 
velocity. 
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Figure 15: The contours of the scalar streamfunction ip(x,y) and the vector 
field (v x .v y ) for a four- vortices solution. 
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The magnitude of the tangential velocity along a radius, v (r) 

200 | 1 1 1 1 1 1 1 r 
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100 - 



50 - 
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r 

Figure 16: The magnitude of the tangential component v$(x,y), seen along 
a radial line. The central dip is visible but significantly narrower than at 
p — 1. 
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3.2 Relevance of the solutions for the physics of the 
atmosphere 



In general the space variables of the CHM equation are normalized to the 
intrinsic typical length of the model. In this case (atmospheric physics) are 
scaled with p g , the Rossby radius. We note in passing, (especially for plasma 
physicists) that there is a major difference compared with the plasma case. 
In plasma, perturbations with lengths less or comparable with an ion Larmor 
gyroradius /c -1 > pi cannot be described by fluid models. 

In the physics of atmosphere, the wavelengths can be much smaller 

kp g > 1 

At very large kp g the description becomes governed by the Euler equation 
(see pj). 

For the numerical studies we choose 

(•'•''; 1J ) ^ [^min; ^max] ^ [Z/min; 2/max] 

= [-0.5,0.5] x [-0.5,0.5] 

This means that the full domain (the side of the rectangle) is a single unit 
length p g . 

In the following we make few consideration about what we can expect as 
results, in the case of the atmosphere problem. 

As we will notice from numerical solution, the equation produces functions 
with very clear similarity with the typhoon morphology. The characteristic 
aspect is (within the precision of these first integrations) a sharp extremum 
of the vorticity on (0, 0) which means a localised maximum of the tangential 
velocity vq in close proximity of the center. Since 

8 dr 

the maximum at 

r = a 

means 

dve d 2 ^) 
dr dr 2 

The equation is 

d 2 ip 1 # ( 1\ ., w 

77 H r- = -777 smh ip (cosh ip - p) 

dr A r dr \ 2p z J 
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We multiply by r and we make a derivation to r 



d 2 ip d 2 dip d 2 ip 



dr 2 



2p 2 



dr 2 dr dr 2 
{sinh ip (cosh ip — p) + 



+r [cosh 2 ip — p cosh ip + sinh 2 ip\ | 



We calculate this expression and the equation in the point r = a defined as 
the point of the maximum of the tangential velocity. This means 



d 2 jj 
dr 2 , a 

dr 2 \ dr 
ip (r = a) 



a 




d 2 v e 



dr 2 



= —a where a > 



where we have introduced a notation for the value, —a < of the second 
derivative of the tangential velocity at its maximum. For a very qualitative 
estimation, used in predicting shapes of solutions, we will take this as a 
parameter. At the point r = a the equation becomes 



1 



dip 
dr 



1 



2p 



sinh ip (cosh ip — p) 



In the equation derivated at r we replace dip/dr with its value from the above 
equation and also introduce the parameter a. Then we have 



-a 



, {sinh ip (cosh ip Q - p) 

Zp 



2p 2 



sinh ip (cosh ip — p) 



x (2 cosh 2 ip — p cosh ip — l) } 



or 



aa (2p 2 ) 



[2 cosh 2 ipQ — p cosh ipQ — l) 



sinh ipQ (cosh ip — p) 2p 2 

This equation may serve to make some estimates if additional informations 
(or simply hints from experiments) are available. This is illustrated below. 
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Consider the case p — 1 



sinh-^o (cosh-^o — 1) 



2aa 




A short and dirty approximation should start by using the suggestion 
from results of lucky simulations, where ipo is few units, and a is of the order 
0.1 on a domain of length 1 in both x and y. The second derivative of the 
tangential velocity must be high, shown by the plots of v$. This means that 
it may exist a difference of magnitude of the terms, with the second term in 
the right hand side appearing less important. Therefore we try 



In addition, we can suppose that the exponentials of negative argument are 
less important than those of positive argument, and simplify to 



2aa ~ sinh ■?/>() (cosh-^ — 1) 



exp (2^ ) ~ 8aa 



or 



■00 ~ ~ In (aa) + 1 



For an order of magnitude we may take 



a ~ — 

a 6 



and then 




In a + 1 



We obtain 



V> - 1 - « In ipo 



— In a 



1 



exp — 1 — - In ipo 
^=exp (4> - 1) 



a 




or 
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We can see that the results are consistent, since if we take from numerical 
solution 

~ 3 



we obtain from the estimation 



which is not far from 



a ~ 0.032 



We must remember that the domain of integration is of length 1 and the fact 
that "the radius of maximum wind" is so small, a ~ 0.04 , means that high 
accuracy is needed to describe correctly what happens close to the center. 
This is due to the other constraint, that the solution streamfunction ip (r) 
needs sufficient space to go to the constant value at "infinity" (large r) . Any 
restriction of the domain of integration which would be aimed to the better 
description of the central region would require boundary conditions that are 
unknown. 

There is another benefit from these very rough estimations. We can use 
them to determine the spatial domain that would be adequate for the search 
of the solution, for particular physical situations. 

In order to use this rough estimation we must introduce physical units. 
In the following all quantities with physical dimensions have an superscript 
phy. 

In atmosphere the distances are measured in p g 

a phy 



and the streamfunction is normalised with 

^phy 



1> 



Pl(f) 

where (/) is the Coriolis parameter. This means 




exp 



r hV 
Pi (f) 



The physical parameters are (taken from |15j ) 
The depth of the atmosphere 



H = 8x 10 3 (m 
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The Coriolis parameter 



(/) = 1.6 x 10- 4 (s- 1 ) 

From these parameters it results 

The Rossby radius (the unit of space) 

_ (gHfl 
99 ~ (f) 

= 2 x 10 6 (m) 
The unit for the streamfunction is 

p\ (f) = 6.4 x 10 8 {m 2 /s) 

The unit for vorticity 

(f) = 1.6 x 10- 4 (s- 1 ) 

For example, using these parameters, it results that we have integrated 
on a spatial domain of length L (in other words: we have imposed that the 
streamfunction becomes equal to ip h on the boundaries of a square with 
side length L) 

L = ^min 1 

L phy = 1 x p g ~ 2 x 10 6 (m) = 2000 (km) 

and the diameter d of the eye of the typhoon results 

d = 2 x a = 0.08 
d phy ~ 0.08p 3 = 128 (km) 

In the Ref. |24j it is reproduced a plot of an observation made on the profile 
of the vorticity, in Fig. la. The plot indicates a maximum value of about 
250 x 10~ 4 (s^ 1 ). The vorticity we obtain is larger (of the order of 1000 x 
10 -4 ). This shows that the absence of the third dimension in our model 
and of the viscous effects have a serious influence on the physical quantities. 
They should be somehow accounted for by renormalizing the two-dimensional 
model at the initial stage. For example, in the case of the plasma vortex, a 
change of the space scale results from the presence of a translational motion 
combined with the density gradient. This remains to be studied. 
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4 Summary 



We again underline that this equation is very difficult to solve, although it 
requires reasonable computer resources. The main problem is the complexity 
of the space of solutions and the need to explore carefully much of this space 
in order to establish the basins of attraction. We are not able at this moment 
to connect in some practically useful way the sharp transitions between the 
attractors with the stability of the solutions. 

It seems that the solution where the streamfunction ip (x, y) is approxi- 
mately radially symmetric, strongly peaked in origin, is a significant attrac- 
tor, at the level of this very sensitive equation. It presents the particularity 
that the vorticity is practically zero for almost all spatial domain and is 
strongly localised, almost singular, close to the maximum. The aspect of 
this solution is very similar to the two-dimensional image of a typhoon. 

We have several arguments in favor of the conclusion that our equation 
(fTj) may represent the hydrodynamic part of the atmospheric vortex. We 
mention some of them. 

1. The profile of the magntitude of the tangential velocity, as represented 
in Fig. 2 of Ref. [22] is very similar to our FiglT2l This is also confirmed 
by similarity with the Fig. la from Ref( |23j ) : 

2. The profile of the vorticity uj shown in our Fig0]is very similar to Fig. la 
from Ref. [23; 

3. We note that in a series of reported numerical simulations, the ten- 
dency of the fields is to evolve toward profiles that are very close to 
those shown in our figures El [U and H21 For example, the Fig. 7a and b 
of Ref. [21] show the evolution of the azimuthal mean of the vorticity 
and tangential velocity from initial profiles which correspond to a nar- 
row ring of vorticity to profiles that show clear ressemblance with our 
figures 0] and ^] or El The same striking evolution to profiles similar to 
ours appears in Figs. 7 a and b of the same Reference. We have investi- 
gated whether a radially annular profile of vorticity can be a solution of 
our equation ([TJ. The result is negative, which may explain why such 
an initial profile evolves to either a set of vortices (vortex-crystal) or 
to a centrally peaked structure as in FiglHl 

4. The four vortices represented in Figure 4a of the Ref. [21] as the late 
stage of the evolution obtained from numerical simulation of vorticity, 
is clearly similar to our figure [T3] 
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5. We obtain a good consistency between our quantitative results for an 
atmospheric vortex (using most elementary input information) and the 
values measured or obtained in numerical simulations, at least for some 
of the quantities. 

A large database on typhoons can be found in [28 . The similarity is 
striking and it suggests that further work with this equation is worth to be 
done. 

The numerical simulations we have taken as a comparison are very com- 
plex. In general, the physics of the typhoons is very complex and includes 
hydrodynamics and thermic aspects, with many additional elements: pre- 
cipitation, viscosity, etc. In no way we do not claim that this equation 
can represent this complexity. It appears however useful as a description of 
the regimes where the hydrodynamical processes are dominating and have 
reached stationarity. 
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5 Appendix A : The structure of a radial so- 
lution near r = and r = oo 

The equation we discuss is 



Other members of the family of equations (parametrized by solutions of the 
2D Laplace equation) will be examined separately. Their importance stems 
from the fact that they can provide, in principle, azimuthal trigonometric 
variation, as for example the Larichev-Reznik modon. 

5.1 The behavior near r = 

Close to the origin, in a purely radial form, it is 



Ai/j H sinh ijj (cosh ip 



p)=0 




H — h I — - I sinh if) (cosh ip — p) = 



r dr 
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where r is measured in p s . 

We take an expansion with only even powers of r close to the origin 

ip ~ a + a 2 r ' 2 + + 06^ 6 + ••• 



Then, for small r; 



^ = 2a 2 r + 4a 4 r 3 + 6a 6 r 5 ... 
ar 

1 dip 2 4 

— — = 2a 2 + 4a 4 r + 6a6?" ... 
r ar 

= 2a 2 + 12a 4 r 2 + 30a 6 r 4 + ... 

ar^ 

sinh (a + a 2 r 2 + Q4 r4 + •••) 
= sinh a 

+ (a 2 r 2 + a 4 r 4 + ...) coshao 

+ \ ( a 2 r4 + •••) sinha + ... 

cosh (ao + a 2 r 2 + a 4 r 4 + ...) 
= cosh a 

+ (a 2 r 2 + a^r 4 + ...) sinh ao 

+- (air 4 + ...) cosha + ••• 



Introducing the notations 



U = a 2 r + a 4 r + ... 
1/ = i(a 2 r 4 + ...) 

sinh ^ (cosh ip — p) 
= (sinh ao + £/ cosh a + sinh ao) 

x (cosh clq — p + U sinh a + V cosh a ) 
= sinh ao (cosh a — p) 

+U (sinh 2 a + cosh 2 ao — p cosh ao) 

+V (2 cosh a sinh a — p sinh a ) 

+£/ 2 (sinh ao cosh ao) 

+V 2 (sinh ao cosh ao) 

+UV (cosh 2 a + sinh 2 a ) 

+... 
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We collect the various degrees of r a 



<?o + (br 2 + W A + ■■■ 
q = sinh a (cosh a — p) 
<& = 0,2 (sinh 2 ao + cosh 2 a — p cosh do) 

g 4 = a 4 (sinh 2 a + cosh 2 a — p cosh a ) 
1 

+ -a 2 (2 cosh a sinh a — p sinh a ) 
+a 2 (sinh a cosh a ) 



Returning to the differential operator 

dr 2 r dr 
= 2a 2 + f2a 4 r 2 + 30a 6 r 4 + ... 

+2a 2 + 4a 4 r 2 + 6a 6 r 4 ... 

= 4a 2 + f6a 4 r 2 + 36a 6 r 4 + ... 

We now identify the expressions corresponding to the same degrees of r, 

4a 2 + 16a 4 r 2 + 36a 6 r 4 + ... 



+ (^) (* + ^ 2 + ^ 4 + ...) 



= 



with the equalities 



4 ° 2 + ) * = ° 

1 6«4+(^]g 2 = o 

36fl 6 + (^2 ) <?4 = 



4a 2 + ( — ) sinh a (cosh ao — p) = 



The equations from which we derive the coefficients of the expansion become 

2p 2 

16a 4 + (^^j a 2 (sffih 2 °o + cosh 2 a — p cosh a ) = 
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36a 6 + (^^T J [ a4 (si 11 ^ 2 «o + cosh 2 ao — J9cosha ) 

+-a2 (2 cosh a sinh a — p sinh a ) 

+0 2 , (sinh a cosh a )] 
= 

We see that if we take 

a = 

then this will vanish all the other coefficients 

a 2 = 
a 4 = 
a 6 = 0,... 




Figure 17: Coefficient a 2 for p = 1. 

Consider the value of the constant 

p — 1 

and we choose the main coefficient of the expansion close to r = to be 

a = 1 

Then 

a 2 = -0.0798 
a A = 0.0055 
a 6 = -0.000439 
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Figure 19: 



Coefficient a 6 for p = 



1 



But the coefficients, as shown in the Figures, are very rapidly growing in 
absolute value. 

We conclude that any attempt to identify the solution starting from few 
terms expansion around r = will be imprecise. 



5.2 The behavior at infinity 

At r — > oo we expect that the function approaches zero in the case where 
p = 1 or approaches one of the roots of the equation 

cosh if) — p = (9) 

for p > 1. The case where if) — > will be treated below. We note, for the 



case p > 1 that the solutions of the Eq.fjHJ are 

In (p + \lp 2 - lj 
In (p — \lp 2 — \ 



5.2.1 The case p=l 

This requires that if) — > at r 
Change the variable 



oo. 



dr 



dr 2 



dx d 
dr dx 

d f d 
dr \ dr 

-x 2 



1 

X 

1 d 

r 2 dx 



-x 



d_ 

dx 



d 



—x 



-X 



d 

-2x— 

dx 



dx dx z 



dx 

x 2 —) 
dx 2 J 



d 

dx 



The function 



— - sinh if) (cosh if) — p) 
2p 2 



if; -> 



2p 2 
1 — p 
~2p T 
1 

f V 



i 

2 



1 — p 



if) 2 



ip 3 + ... 
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For 



p = 1 

— - ) sinh ip (cosh ib — p) — > — -ip 3 
2p 2 I ' 4 



Then the equation becomes 



dx dx A 



1 - p , 1/1 1 - p\ , 



2p 2 r 2 p 2 V 2 6 





This can be approximated at 

x -> 



or 



For 



-a; 2 ^ = + /30 3 
ax 

a — = dr 



aip + (3tp 3 x 2 \x 

p = 1 
a = 

""4 



dip 

'2 



then 





V r 

We note however that in this case the vorticity is 

uj = Aip 

We would like to have a vanishing vorticity at infinity with a faster decay. 

The above calculations seem to suggest that for purely radial structure 
we need to consider the differential equation which is derived for a different 
choice of the Laplacean equation, as it is explained in the main text. 
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5.2.2 The case p > 1 



One possibility, for 



p > 1 



■0 ~ exp (— \a\ r) 




This gives 



Aip 




exp (— |a| r) 



+ a 2 exp (— |a| r) 



r 



with a fast decay. This situation is worth to be examined numerically. 

6 Appendix B : various forms of the initial 
conditions 

6.1 The ring- type 

The initial form of the function has the form 



and we take the maximum to be placed at 

r = a 

which is considered to approximate the center line of the ring. The equation 
becomes 



ipo = A exp (— sr 2 ) [l — /texp (— gr 4 )] 



We look for the maximum 



dipo 



(— 2sr) exp (— sr 2 ) [l — Kexp (— gr 4 )] 
+ exp (— sr 2 ) (4gr 3 ) /texp (— gr 4 ) 



dr 



= 



(2sk + Anqa 2 } exp (— ga 4 ) = 2s 
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The other condition is that the maximum of the function ip at r = a equals 
a prescribed value, 



(r = a) = ip c 
Aexp (— sa 2 ) [l — /texp (— ga 4 )] = ^ 

The initial condition is introduced in the following way. We take q, a, k 
and ip c as input parameters and determine the other two, s and A from the 
equations 

2a 2 q 

s = — 

Kexp (— ga 4 ) — 1 



exp (— sa 2 ) [1 — /texp (— ga 4 )] 
Now the initial function will be 



^initial (r) = 1p + Vf' 2) 



= i>r + 

+Aexp (— sr 2 ) [l — /texp (— gr 4 )] 

ie. the function just determined is placed on the constant background of 
the value at the boundary, calculated form the condition that the vorticity 
is zero at infinity. 

This class of initial functions is characterised by an annular shape, with 
exponential decay for r — > oo, with a minimum in the region around r = of 
depth that can be fixed by varying k. For k — 1 the function is zero on the 
symmetry axis and rises slowly (due to r 4 ) toward the maximum at r = a. 

In order to narrow the space of parameters we require the approximative 
equality between the vorticity amplitude at the ring with the nonlinear term 

2 , 

u j 2 ^ c 

~ ~2p2 Sinh + COsh + ~ p 

(Here 5 is the width of the ring shape). These two quantities are compared 
in graphical plot for a range of values of the parameter ip c , using a Matlab 
script. This is far from an exact procedure but helps to generate reasonable 
ranges for the input parameters. 

The conclusion after many trials using this procedure and its initial func- 
tion forms can be described as follows. 
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In most of the cases the central region is corrected and shifted to a maxi- 
mum. In the cases p — 1 the central region which is started with a deppressed 
level is rised and a strong peaked form is generated, as in the cases where 
the initialization consists of a maximum on center (for example a Gaussian 
form). For p > 1 the run evolves in some cases to the formation of separate 
maxima placed symmetrically on a ring, having sharp maxima. The central 
region is decreased in amplitude to a somehow flat region. The region outside 
the ring is evolving to a state which corresponds with very good precision, 
to 

uj ~ 

on the rest of the domain to the periphery. 

6.2 Flat central region for ip (r) 

We take the central region 

< r < r flat 

with a fixed, constant value 

i> (r) = i) c 

where ip c is one of the roots of the equation cosh-?/' — p = 0. At the edge we 
take another fixed value, 

i> = ipb 

with ipb the other, smaller root of the equation. 
In between, we take 

■0 (r) = -01 — A In (r) 

A _ ^c~A 
In (r f iat/r c ) 
4>i = i> c + Aha. (r f iat) 

The value r c is 

rfiat <r c < 0.5 

represents the value where where we stop the decay of the function with 
logarithm profile and put ip = ipb- This is 

rfiat = 0.1 

r c = 0.35 •••0.45 

The parameter p = 1.3. 
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The result of these calculations is as follows. For small mesh, the evolution 
is clearly toward the suppression of the smoothly decaying part, letting a sort 
of cylinder in the center, with radius rfi at , with the high value equal to ip c 
and the rest seems to go progressively to ip — ipj,. The vorticity is singular, 
around r = r// at . The vorticity is positive and negative, with high values, 
singular in a narrow ring. 

For this cylindrical-rod profile of the streamfunction ip (r), the velocity is 
very localised, as a very narrow ring, all its values are positive. The velocity 
grows from zero, keeps always the same direction on 9 and then decays to 
zero value, after the width of the ring. The vorticity is also sharply limitted 
here, but it has positive and negative values on interior half of the ring and 
respectively on the exterior half of the ring. 

The same shrinking to the cylindrical column happens when we take the 
maximum of tp (in the central flat region) as 

ip (r) = 

which is the other possibility that the equation is verified for constant value 
of ip. 
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